######################################################
#### Replication Material for II Joo & Mukherjee #####
###     Manuscript ID: GINI-2019-1752.R2 (new)   #####
####                   GINI-2019-1752.R1 (old)   #####  
####             Main Paper Analyses             #####
####       Software Used: R Version 3.6.0        #####
######################################################
rm(list=ls())


### use for creating a log file ###

# wd <- setwd("Your Directory/Joo_Mukherjee_Replication")
# (wd <- getwd())

# sink("joo_mukherjee_main_R_output.log")

### set directory ###

wd <- setwd("Your Directory/Joo_Mukherjee_Replication")
(wd <- getwd())


### load libraries ###

library(foreign)
library(stargazer)
library(ggplot2)
library(mvtnorm)
library(ggthemes)
library(reshape2)
library(margins)
library(survival)
library(survminer)



jm <- read.dta("C:/Users/happy/Dropbox/Penn State/Dissertation/splits/Joo_Mukherjee_Replication/II_Joo_Mukherjee_12.dta")

################
### Figure 1 ###
################

coef <- c(0.4455619,
          0.0157448,
          -0.3074662,
          -0.0343053,
          0.264737,
          0.2530055,
          0.2918376,
          -0.1933579,
          0.0040991,
          -0.2345839,
          -0.0237582,
          -0.0040464,
          -27.90933,
          2.099829,
          -0.0496096,
          -0.9857114)

vcov <- as.matrix(read.csv("vcov.csv"))

lag <- 0
rebst <- 1
nego <- 0
cease <- 0
mob <- 1
lbd <- median(jm$logbd, na.rm=T)
lbd2 <- median(jm$logbd2, na.rm=T)
ind <- 0
terri <- 0
t <- median(jm$t1, na.rm=T)
t2 <- median(jm$t2, na.rm=T)
t3 <- median(jm$t3, na.rm=T)
censt.sim <- c(1,2,3)
mtime <- median(jm$dur_exist, na.rm=T)
sdtime <- sd(jm$dur_exist, na.rm=T)

set.seed(88110905)
B <- 3000
out <- matrix(NA, ncol=3, nrow=B)
ext <- matrix(NA, ncol=6, nrow=B)

for(i in 1:length(censt.sim)){ 
  censt <- censt.sim[i] 	
  pr.bs <- matrix(NA,B,1)
  for(b in 1:B){ 
    newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
    bs.sum <- newcoef[1]*lag + newcoef[2]*censt.sim[i]*(mtime) + newcoef[3]*censt.sim[i] + newcoef[4]*(mtime) + newcoef[5]*nego + 
      newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*ind + newcoef[11]*lbd + newcoef[12]*lbd2 +
      newcoef[13]*t + newcoef[14]*t2 + newcoef[15]*t3 + newcoef[16]
    
    
    bs.sum2 <- newcoef[1]*lag + newcoef[2]*censt.sim[i]*(mtime+sdtime) + newcoef[3]*censt.sim[i] + newcoef[4]*(mtime+sdtime) + newcoef[5]*nego + 
      newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*ind +  newcoef[11]*lbd + newcoef[12]*lbd2 + 
      newcoef[13]*t + newcoef[14]*t2 + newcoef[15]*t3 + newcoef[16]
    
    out[b,i] <- pnorm(bs.sum2) - pnorm(bs.sum)
    
    if(i==1){
      ext[b,1] <-pnorm(bs.sum2)
      ext[b,2] <- pnorm(bs.sum)
    } else{
      if(i==2){
        ext[b,3] <-pnorm(bs.sum2)
        ext[b,4] <- pnorm(bs.sum)
      } else{
        if(i==3){
          ext[b,5] <-pnorm(bs.sum2)
          ext[b,6] <- pnorm(bs.sum)
        }
      }
    }
    
  }
  rm(bs.sum)
  rm(bs.sum2)
}


low <- as.matrix(sort(out[,1]))
moder <- as.matrix(sort(out[,2]))
high <- as.matrix(sort(out[,3]))

sort.out <- cbind(low, moder, high)
conf95 <- sort.out[(B*0.05+1):(B*0.95),]


pdf(paste(wd, "/Figure_files/Figure1.pdf", sep=""),width=8,height=8,paper='special') 

boxplot(conf95, boxwex=0.5, names=c("Low", "Moderate", "High"), xlab="Central Command", 
        ylab="Change in Pr(Split)", cex.lab=1.5, cex.axis=1.5, cex=1.5)
abline(h=0, lty=2, col='red')

dev.off()


#####################
#### Figure 2a-c ####
#####################

##### set up #####

coef <- c(0.3732367,
          0.0169961,
          -0.355772,
          -0.0414456,
          0.2772929,
          0.4089621,
          0.3525784,
          -0.1484045,
          -0.0332967,
          -37.63297,
          3.684742,
          -0.0841062,
          -1.430015)

vcov <- as.matrix(read.csv("vcov_fig2.csv"))


lag <- 0
rebst <- 1
nego <- 0
cease <- 0
mob <- 1

ind <- 0
terri <- 0
t <- median(jm$t1, na.rm=T)
t2 <- median(jm$t2, na.rm=T)
t3 <- median(jm$t3, na.rm=T)
censt.sim <- c(1,2,3)
mtime <- median(jm$dur_exist, na.rm=T)
sdtime <- sd(jm$dur_exist, na.rm=T)
time.sim <- seq(min(jm$dur_exist), max(jm$dur_exist), by=0.5)



#############################
#### Figure 2a: High  #######
#############################
B <- 1000

set.seed(905)
out1 <- matrix(NA, ncol=3, nrow=length(time.sim))
for(i in 1:length(time.sim)){ 
  time <- time.sim[i] 	
  pr.bs <- matrix(NA,B,1)
  for(b in 1:B){ 
    newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
    bs.sum <- newcoef[1]*lag + newcoef[2]*3*time + newcoef[3]*3 + newcoef[4]*time + newcoef[5]*nego + 
      newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*t + newcoef[11]*t2 + newcoef[12]*t3 + newcoef[13]
    
    
    pr.bs[b,] <- pnorm(bs.sum) 
  }
  out1[i,] <- c(mean(pr.bs, na.rm=TRUE), quantile(pr.bs, c(0.025, 0.975), na.rm=TRUE)) 
  rm(bs.sum)
}

pdf(paste(wd, "/Figure_files/Figure2a.pdf", sep=""),width=10,height=8,paper='special') 

par(mar=c(5,5,1,1))
plot(time.sim, out1[,1], type="n", xlab="Rebel Age", ylab="Pr(Split)", 
     ylim=c(-0.0005,0.3),
     cex.lab=1.8, cex.axis=1.5)
lines(lowess(time.sim, out1[,1]), lty=1, lwd=2)
lines(lowess(time.sim, out1[,2]), lty=2, lwd=2)
lines(lowess(time.sim, out1[,3]), lty=2, lwd=2)

dev.off()


##############################
#### Figure 2b: Moderate  ####
##############################

B <- 1000

set.seed(905)
out3 <- matrix(NA, ncol=3, nrow=length(time.sim))
for(i in 1:length(time.sim)){ 
  time <- time.sim[i] 	
  pr.bs <- matrix(NA,B,1)
  for(b in 1:B){ 
    newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
    bs.sum <- newcoef[1]*lag + newcoef[2]*2*time + newcoef[3]*2 + newcoef[4]*time + newcoef[5]*nego + 
      newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*t + newcoef[11]*t2 + newcoef[12]*t3 + newcoef[13]
    
    pr.bs[b,] <- pnorm(bs.sum)
  }
  out3[i,] <- c(mean(pr.bs, na.rm=TRUE), quantile(pr.bs, c(0.025, 0.975), na.rm=TRUE)) 
  rm(bs.sum)
}

pdf(paste(wd, "/Figure_files/Figure2b.pdf", sep=""),width=10,height=8,paper='special') 
par(mar=c(5,5,1,1))
plot(time.sim, out3[,1], type="n", xlab="Rebel Age", ylab="Pr(Splits)", 
     ylim=c(0.0005,0.15),
     cex.lab=1.8, cex.axis=1.5)
lines(lowess(time.sim, out3[,1]), lty=1, lwd=2)
lines(lowess(time.sim, out3[,2]), lty=2, lwd=2)
lines(lowess(time.sim, out3[,3]), lty=2, lwd=2)

dev.off()


########################
#### Figure 2c: Low ####
########################
B <- 1000

set.seed(905)
out2 <- matrix(NA, ncol=3, nrow=length(time.sim))
for(i in 1:length(time.sim)){ 
  time <- time.sim[i] 	
  pr.bs <- matrix(NA,B,1)
  for(b in 1:B){ 
    newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
    bs.sum <- newcoef[1]*lag + newcoef[2]*1*time + newcoef[3]*1 + newcoef[4]*time + newcoef[5]*nego + 
      newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob +newcoef[10]*t + newcoef[11]*t2 + newcoef[12]*t3 + newcoef[13]
    
    pr.bs[b,] <- pnorm(bs.sum) 
  }
  out2[i,] <- c(mean(pr.bs, na.rm=TRUE), quantile(pr.bs, c(0.025, 0.975), na.rm=TRUE)) 
  rm(bs.sum)
}

pdf(paste(wd, "/Figure_files/Figure2c.pdf", sep=""),width=10,height=8,paper='special') 
par(mar=c(5,5,1,1))
plot(time.sim, out2[,1], type="n", xlab="Rebel Age", ylab="Pr(Split)", 
     ylim=c(-0.0005,0.15),
     cex.lab=1.8, cex.axis=1.5)
lines(lowess(time.sim, out2[,1]), lty=1, lwd=2)
lines(lowess(time.sim, out2[,2]), lty=2, lwd=2)
lines(lowess(time.sim, out2[,3]), lty=2, lwd=2)
dev.off()


##############################################
#### Figure 3: Moderate - High Time Plot #####
##############################################


coef <- c(0.4455619,
          0.0157448,
          -0.3074662,
          -0.0343053,
          0.264737,
          0.2530055,
          0.2918376,
          -0.1933579,
          0.0040991,
          -0.2345839,
          -0.0237582,
          -0.0040464,
          -27.90933,
          2.099829,
          -0.0496096,
          -0.9857114)

vcov <- as.matrix(read.csv("vcov.csv"))


lag <- 0
rebst <- 1
nego <- 0
cease <- 0
mob <- 1
lbd <- median(jm$logbd, na.rm=T)
lbd2 <- median(jm$logbd2, na.rm=T)
ind <- 0
terri <- 0
t <- median(jm$t1, na.rm=T)
t2 <- median(jm$t2, na.rm=T)
t3 <- median(jm$t3, na.rm=T)
censt.sim <- c(1,2,3)
mtime <- median(jm$dur_exist, na.rm=T)
sdtime <- sd(jm$dur_exist, na.rm=T)

time.sim <- seq(min(jm$dur_exist), max(jm$dur_exist), by=0.5)

##### number of bootstraps ##### 
B <- 1000

set.seed(905)
out2 <- matrix(NA, ncol=3, nrow=length(time.sim))
for(i in 1:length(time.sim)){ 
  time <- time.sim[i] 	
  pr.bs <- matrix(NA,B,1)
  for(b in 1:B){ 
    newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
    bs.sum <- newcoef[1]*lag + newcoef[2]*3*time + newcoef[3]*3 + newcoef[4]*time + newcoef[5]*nego + 
      newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*ind + newcoef[11]*lbd + newcoef[12]*lbd2 +
      newcoef[13]*t + newcoef[14]*t2 + newcoef[15]*t3 + newcoef[16]
    bs.sum2 <- newcoef[1]*lag + newcoef[2]*2*time + newcoef[3]*2 + newcoef[4]*time + newcoef[5]*nego + 
      newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*ind + newcoef[11]*lbd + newcoef[12]*lbd2 +
      newcoef[13]*t + newcoef[14]*t2 + newcoef[15]*t3 + newcoef[16]
    pr.bs[b,] <- pnorm(bs.sum) - pnorm(bs.sum2)
  }
  out2[i,] <- c(mean(pr.bs, na.rm=TRUE), quantile(pr.bs, c(0.05, 0.95), na.rm=TRUE)) 
  rm(bs.sum)
}

pdf(paste(wd, "/Figure_files/Figure3.pdf", sep=""),width=12,height=8,paper='special') 
par(mar=c(5,5,1,1))
plot(time.sim, out2[,1], type="n", xlab="Rebel Age (years)", ylab="Change in Pr(Splits)", 
     ylim=c(-0.03,0.15),
     cex.lab=1.8, cex.axis=1.5)
lines(lowess(time.sim, out2[,1]), lty=1, lwd=2)
lines(lowess(time.sim, out2[,2]), lty=2, lwd=2)
lines(lowess(time.sim, out2[,3]), lty=2, lwd=2)
abline(h=0, lty=2, col='red', lwd=1)

dev.off()

# end of file 

# sink() # end log